--- title: "Manejo de imágenes" date: 2020-02-10T17:34:30-04:00 categories: - blog tags: - post - r - spatial - raster layout: splash toc: true ---
Para descargarla desde mi googleDrive vamos a utilizar los siguientes comandos
## Loading required package: sp
library(curl)
#Ubicación para guardar la imagen descargada
location <- "D:/Descargas/imagen1.tif"
#Copiar el link para compartir desde google drive
GD_share_URL = "https://drive.google.com/open?id=17PyQnEIICpjNtPP3v59gLEaoY3zzzEcR"
# Reemplazar "open?" con "us?export=download&"
dwnld_URL <- gsub("open\\?", "uc\\?export=download\\&", GD_share_URL )
dl <- curl_download(dwnld_URL,
destfile = location)
#Leer imagen
imagen1 <- raster(location)Recordar que:
Cargar paquete raster
Escribir completa la ubicación del archivo
Incluir extension
Para ubicaciones de archivos “/” en lugar de “’'”
A continuación, se puede graficar como cualquier gráfica en R
Sin embargo, aquí nada más se cargó una banda. Para cargar una imagen con más de una banda hay que utilizar otra función.
También se puede utilizar la función brick
Regresamos a la imagen importada con stack
Ver las características del objeto ¿Qué es cada cosa?
## class : RasterStack
## dimensions : 81, 103, 8343, 10 (nrow, ncol, ncell, nlayers)
## resolution : 0.0001796631, 0.0001796631 (x, y)
## extent : -101.2353, -101.2168, 19.64364, 19.65819 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
## names : B2, B3, B4, B5, B6, B7, B8, B8A, B11, B12
## min values : 0, 0, 0, 0, 0, 0, 0, 0, 0, 0
## max values : 65535, 65535, 65535, 65535, 65535, 65535, 65535, 65535, 65535, 65535
## [1] 10
## CRS arguments:
## +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
Más características
## [1] 0.0001796631 0.0001796631
## [1] 8343
## [1] 81 103 10
Más características
## [1] "B2" "B3" "B4" "B5" "B6" "B7" "B8" "B8A" "B11" "B12"
## class : Extent
## xmin : -101.2353
## xmax : -101.2168
## ymin : 19.64364
## ymax : 19.65819
¿Qué es esto?
## class : RasterLayer
## band : 1 (of 10 bands)
## dimensions : 81, 103, 8343 (nrow, ncol, ncell)
## resolution : 0.0001796631, 0.0001796631 (x, y)
## extent : -101.2353, -101.2168, 19.64364, 19.65819 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
## source : D:/Descargas/imagen1.tif
## names : B2
## values : 0, 65535 (min, max)
¿y esto?
## [1] "INT2U" "INT2U" "INT2U" "INT2U" "INT2U" "INT2U" "INT2U" "INT2U" "INT2U"
## [10] "INT2U"
Imágenes mutiespectrales con 13 bandas de distinta resolución espacial. Por eso, en este stack sólo se tienen las imágenes
También se puede ver la imagen con otros paquetes
## Loading required package: lattice
## Loading required package: latticeExtra
##
## Attaching package: 'ggplot2'
## The following object is masked from 'package:latticeExtra':
##
## layer
Otra forma de mostrar imagenes
## Linking to GEOS 3.6.1, GDAL 2.2.3, PROJ 4.9.3
## Warning: Raster values found that are outside the range [0, 2100]
Se puede realizar de dos maneras
Primero cargar la imagen y calcular el NDVI
library(raster)
#Ubicación para guardar la imagen descargada
location_dem <- "D:/Descargas/img_DEM.tif"
#Copiar el link para compartir desde google drive
GD_share_URL = "https://drive.google.com/open?id=1qOJxEaAFgYMr8PQNuIKJgcTYq_SyvhO0"
# Reemplazar "open?" con "us?export=download&"
dwnld_URL <- gsub("open\\?", "uc\\?export=download\\&", GD_share_URL )
dl <- curl_download(dwnld_URL,
destfile = location_dem)
dem<-raster(location_dem)
plot(dem)## CRS arguments:
## +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
## [1] 0.0002694946 0.0002694946
## [1] 55 69 1
Otras formas de realizar lo mismo
EVI <- 2.5 * (imagen1$NIR - imagen1$R) / ((imagen1$NIR + 6 * imagen1$R - 7.5 * imagen1$B + 1))
plot(EVI)No se alcanza a distinguir nada. Se ven valores de -40 a 20. Probablemente hicimos algo mal. A veces es mejor pasar de la escala 0 - 10000 a una de 0 - 1 antes de hacer los cálculos de algunos de estos índices. Entonces dividimos cada banda entre 10000.
EVI <- 2.5 * (imagen1$NIR / 10000 - imagen1$R/ 10000) / ((imagen1$NIR/ 10000 + (6 * imagen1$R/ 10000) - (7.5 * imagen1$B/ 10000) + 1))
plot(EVI)Para calcular este índice requerimos otras bandas: B8 y B11, así que volvemos a cargar la imagen.
imagen1<-stack(location)
imagen1<-imagen1/10000
NDWI <- (imagen1[[8]] - imagen1[[9]]) / (imagen1[[8]] + imagen1[[9]])
plot(NDWI)#Tasselled Cap - Brightness
#Ahorita no se va a calcular porque no tenemos la B10
# Brightness <- 0.3037*imagen1[[1]]+0.2793*imagen1[[2]] +0.4743*imagen1[[3]]+0.5585*imagen1[[4]]+ 0.5082*B10+0.1863*imagen1[[10]]
#Tasselled Cap - Greeness
Greeness <- (-0.2848*imagen1[[1]])+(-0.2435*imagen1[[2]])+(-0.5436*imagen1[[3]]) +0.7243*imagen1[[8]]+0.0840*imagen1[[9]]+(-0.1800*imagen1[[10]])
#Tasselled Cap - Wetness
Wetness <- 0.1509*imagen1[[1]]+0.1973*imagen1[[2]]+0.3279*imagen1[[3]] +0.3406*imagen1[[8]]+(-0.7112*imagen1[[9]])+(-0.4572*imagen1[[10]])
plot(Greeness)Este método se basa en la matriz de co-ocurrencias de tonos de gris (GLCM). Calcula diferentes métricas en ciertas direcciones. Es un método de “moving window”
Primero cargar la imagen y calcular el NDVI
library(raster)
library(glcm)
imagen1<-stack(location)
imagen1<-subset(imagen1, c(1:3,8))
names(imagen1)<-c("B","G","R","NIR")
NDVI <- (imagen1$NIR - imagen1$R) / (imagen1$NIR + imagen1$R)Aquí podemos indicar el tamaño de la ventana que vamos a utilizar para calcular la textura. Para que sólo un pixel obtenga el valor de la textura de la ventana hay que utilizar un tamaño non. El parámetro shift va decir cuántos pixeles se va a mover en X y en Y. En este caso el shift indica que se va a comparar cada pixel dentro de la ventana con el pixel que esté a una distancia de un pixel a la derecha y uno arriba.
glcm_NE <- glcm(NDVI,
window = c(9,9),
shift = c(1,1),
statistics = c("mean", "variance", "homogeneity", "contrast",
"dissimilarity", "entropy", "second_moment", "correlation")
)
plot(glcm_NE)Debido a que se puede calcular la textura en 4 direcciones, para obtener una medida sin efecto de dirección se ponen las cuatro posibilidades: (0,1); (1,1); (1,0); (1,-1)
glcm_Alldir<- glcm(NDVI,
window = c(9,9),
shift=list(c(0,1), c(1,1), c(1,0), c(1,-1)),
statistics = c("mean", "variance", "homogeneity", "contrast",
"dissimilarity", "entropy", "second_moment")
)
plot(glcm_Alldir)Este método primero realiza una transformada de Fourier de la imagen y luego realiza una ordenación de estos datos. Por eso se llama Ordenación de la Transformada de Fourier (FOTO). La idea de este método es que caracteriza la textura a partir de la frecuencia de las ondas dominantes (r-spectrum). De tal manera, si una ventana presenta una textura dominante a una distancia de varios pixeles, estará caracterizado por ondas de frecuencia corta (pocos ciclos por km). En cambio, si una ventana presenta una textura dominante a una distancia de pocos pixeles, estará caracterizado por ondas de frecuencia largas (muchos ciclos por km). Este método se calcula por áreas de la imagen que correponden al tamaño de la ventana.
Este método sólo permite calcularse para ventanas cuadradas. Por ello, sólo se indica el lado de la ventana cuadrada.